/****************************************************************************
*                                                                           *
* scatblock.cpp																*
*                                                                           *
* ALaRI 2008												                *
*                                                                           *
****************************************************************************/

/****************************************************************************
* Includes										                            *
****************************************************************************/
#include "scatblock.h" // - me

/****************************************************************************
* class CTransform - transform helper class									*
****************************************************************************/

/****************************************************************************
* General helpers															*
****************************************************************************/
CTransform::eTransform CTransform::combineTransforms(
	eTransform nTransform1, eTransform nTransform2)
{
	// validate params
	assert(nTransform1 >= eTransform_r0 && nTransform1 < eTransform_count);
	assert(nTransform2 >= eTransform_r0 && nTransform2 < eTransform_count);

	// we could implement some smarter code but we prefer just table
	// NOTE: we know that matrixes are symetric and nice ... but we still ignore this fact!
	static const BYTE s_rgnbyTransforms[eTransform_count][eTransform_count] =
	{
		{ // nTransform1 = eTransform_r0
		eTransform_r0, eTransform_r90, eTransform_r180, eTransform_r270,
		eTransform_fr0, eTransform_fr90, eTransform_fr180, eTransform_fr270
		},
		{ // nTransform1 = eTransform_r90
		eTransform_r90, eTransform_r180, eTransform_r270, eTransform_r0,
		eTransform_fr270, eTransform_fr0, eTransform_fr90, eTransform_fr180
		},
		{ // nTransform1 = eTransform_r180
		eTransform_r180, eTransform_r270, eTransform_r0, eTransform_r90,
		eTransform_fr180, eTransform_fr270, eTransform_fr0, eTransform_fr90
		},
		{ // nTransform1 = eTransform_r270
		eTransform_r270, eTransform_r0, eTransform_r90, eTransform_r180,
		eTransform_fr90, eTransform_fr180, eTransform_fr270, eTransform_fr0
		},
		{ // nTransform1 = eTransform_fr0
		eTransform_fr0, eTransform_fr90, eTransform_fr180, eTransform_fr270,
		eTransform_r0, eTransform_r90, eTransform_r180, eTransform_r270
		},
		{ // nTransform1 = eTransform_fr90
		eTransform_fr90, eTransform_fr180, eTransform_fr270, eTransform_fr0,
		eTransform_r270, eTransform_r0, eTransform_r90, eTransform_r180
		},
		{ // nTransform1 = eTransform_fr180
		eTransform_fr180, eTransform_fr270, eTransform_fr0, eTransform_fr90,
		eTransform_r180, eTransform_r270, eTransform_r0, eTransform_r90
		},
		{ // nTransform1 = eTransform_fr270
		eTransform_fr270, eTransform_fr0, eTransform_fr90, eTransform_fr180,
		eTransform_r90, eTransform_r180, eTransform_r270, eTransform_r0
		}
	};
	// just return from table
	return (eTransform)s_rgnbyTransforms[nTransform1][nTransform2];
}

/****************************************************************************
* CScatBlockL0 - private class for scattering matrix (level 0)				*
****************************************************************************/

/****************************************************************************
* Construction/Destruction													*
****************************************************************************/
CScatBlockL0::CScatBlockL0(const std::complex<float>& fltThetta) :
	m_fltThetta(fltThetta)
{
	// nothing to do
}

/****************************************************************************
* Overrides																	*
****************************************************************************/
int CScatBlockL0::getWidth() const
{
	// only one column
	return 1;
}

int CScatBlockL0::getHeight() const
{
	// only one row
	return 1;
}

const CScatBlock* CScatBlockL0::getTransformed(
	CTransform::eTransform /*nTransform*/, bool /*bCreateIfNotThere*/) const
{
	// we don't need any rotations flips
	return this;
}

std::complex<float>& CScatBlockL0::valueAt(int nRow, int nColumn)
{
	// valid index is only 0,0
	assert(nRow == 0 && nColumn == 0);
	return m_fltThetta;
}

/****************************************************************************
* CScatBlockLP - scat matrix block data for LP matrices						*
****************************************************************************/

/****************************************************************************
* Overrides																	*
****************************************************************************/
int CScatBlockLP::getWidth() const
{
	return m_nWidth; // just getter
}

int CScatBlockLP::getHeight() const
{
	return m_nHeight; // just getter
}

const CScatBlock* CScatBlockLP::getTransformed(
	CTransform::eTransform nTransform, bool bCreateIfNotThere) const
{
	assert(nTransform >= CTransform::eTransform_r0 && nTransform < CTransform::eTransform_count);

	// check null transformation first
	if (nTransform == CTransform::eTransform_r0)
		return this;
	// create if not created
	else if (m_rgnpTransforms[nTransform - 1] == NULL && bCreateIfNotThere)
	{
		CScatBlockLP* pThis = (CScatBlockLP*)this;
		pThis->m_rgnpTransforms[nTransform - 1] = new CScatBlockLPT(pThis, nTransform);
	}

	return m_rgnpTransforms[nTransform - 1];
}

std::complex<float>& CScatBlockLP::valueAt(int nRow, int nColumn)
{
	assert(nRow >= 0 && nRow < m_nHeight);
	assert(nColumn >= 0 && nColumn < m_nWidth);

	// special case for symetric matrixes
	if (m_bSymetric)
	{
		// below diagonal?
		if (nColumn < nRow)
			return m_ppfltData[nColumn][nRow - nColumn];
		else // above diagonal
			return m_ppfltData[nRow][nColumn - nRow];
	}
	else // normal matrix
		return m_ppfltData[nRow][nColumn];
}

/****************************************************************************
* Construction																*
****************************************************************************/
CScatBlockLP::CScatBlockLP(bool bSymetric, int nWidth, int nHeight) :
	m_bSymetric(bSymetric),
	m_nWidth(nWidth),
	m_nHeight(nHeight),
	m_ppfltData(NULL)
{
	// dimentions should be valid
	assert(m_nWidth >= 1 && m_nHeight >= 1);

	// symetric block is always square
	assert(!m_bSymetric || m_nWidth == m_nHeight);

	// reset pointers
	memset(&m_rgnpTransforms[0], 0, sizeof(m_rgnpTransforms));

	// allocate data
	m_ppfltData = allocData(m_nWidth, m_nHeight, m_bSymetric);
}

CScatBlockLP::~CScatBlockLP()
{
	// get rid of data
	freeData(m_ppfltData, m_nHeight, m_bSymetric);

	// delete transforms
	for (int nTransform = 0;
		nTransform < sizeof(m_rgnpTransforms) / sizeof(m_rgnpTransforms[0]); ++nTransform)
	{
		if (m_rgnpTransforms[nTransform] != NULL)
			delete m_rgnpTransforms[nTransform];
	}
}

/****************************************************************************
* Public access																*
****************************************************************************/
void CScatBlockLP::invert()
{
	// we invert only non symetric matrices but square matrices
	assert(!m_bSymetric && m_nWidth == m_nHeight);

	// allocate another copy of data
	std::complex<float> ** ppfltData = allocData(m_nWidth, m_nHeight, false);

	int nRover; // fill with unitary matrix
	for (nRover = 0; nRover < m_nWidth; ++nRover)
		ppfltData[nRover][nRover] = 1.f;

	// calculate using jordan-gauss with pivoting
	for (nRover = 0; nRover < m_nWidth; ++nRover)
	{
		int nRow, nPivotRow; // search for the pivot row starting from current rover
		float fltAbsMax, fltAbs;
		for (nPivotRow = nRover, fltAbsMax = std::abs(m_ppfltData[nRover][nRover]),
			nRow = nRover + 1; nRow < m_nHeight; ++nRow)
		{
			if ((fltAbs = std::abs(m_ppfltData[nRow][nRover])) > fltAbsMax)
			{
				fltAbsMax = fltAbs;
				nPivotRow = nRow;
			}
		}

		// pivot should not be 0.f or the matrix is not invertable
		assert(fltAbsMax != 0.f);

		// swap rows if needed in data matrix and unitary matrix
		if (nPivotRow != nRover)
		{
			std::complex<float> *pfltDataRow = m_ppfltData[nPivotRow];
			m_ppfltData[nPivotRow] = m_ppfltData[nRover];
			m_ppfltData[nRover] = pfltDataRow;

			pfltDataRow = ppfltData[nPivotRow];
			ppfltData[nPivotRow] = ppfltData[nRover];
			ppfltData[nRover] = pfltDataRow;
		}

		// devide pivot row by pivot value in both matrices
		int nColumn; std::complex<float> fltcTemp;
		
		for (nColumn = 0, fltcTemp = m_ppfltData[nRover][nRover]; nColumn <= nRover; ++nColumn) 
			ppfltData[nRover][nColumn] /= fltcTemp;

		for (m_ppfltData[nRover][nRover] = 1.f; nColumn < m_nWidth; ++nColumn)
		{
			m_ppfltData[nRover][nColumn] /= fltcTemp;
			ppfltData[nRover][nColumn] /= fltcTemp;
		}

		// rest are adjusted by pivot column value multiplied by pivot row
		for (nRow = 0; nRow < m_nHeight; ++nRow) 
			if (nRow != nRover)
			{
				// for column prior to rover we need to handle only unitary part
				for (nColumn = 0, fltcTemp = m_ppfltData[nRow][nRover]; nColumn <= nRover; ++nColumn) 
					ppfltData[nRow][nColumn] -= fltcTemp * ppfltData[nRover][nColumn];

				// for all other columns substruct rover column 
				// multiplied by the element at rover column
				for (m_ppfltData[nRow][nRover] = 0.f; nColumn < m_nWidth; ++nColumn) 
				{
					m_ppfltData[nRow][nColumn] -= fltcTemp * m_ppfltData[nRover][nColumn];
					ppfltData[nRow][nColumn] -= fltcTemp * ppfltData[nRover][nColumn];
				}
			}
	}

	// replace data on 'this'
	freeData(m_ppfltData, m_nHeight, m_bSymetric);
	m_ppfltData = ppfltData;
}

/****************************************************************************
* Protected helpers															*
****************************************************************************/
std::complex<float>** CScatBlockLP::allocData(int nWidth, int nHeight, bool bSymetric)
{
	// dimentions should be valid
	assert(nWidth >= 1 && nHeight >= 1);

	// symetric item should be square
	assert(!bSymetric || nWidth == nHeight);

	// number of complex floats needed
	int nSize = bSymetric ? (((nHeight + 1) * nHeight) / 2) : (nHeight * nWidth);
	
	// allocate data for pointers to rows and 'complex float' data
	std::complex<float>** ppfltData = (std::complex<float>**)
		new BYTE[sizeof(std::complex<float>*) * nHeight + nSize * sizeof(std::complex<float>)];

	// construct floats array (this should not throw) - floats come after pointers
	new (&ppfltData[nHeight]) std::complex<float> [nSize];

	int nRow; // now set pointers to rows
	// if symetric - decrement number of elements in each succesive row
	for (nRow = 1, ppfltData[0] = (std::complex<float>*)&ppfltData[nHeight];
		nRow < nHeight; ++nRow, nWidth -= bSymetric ? 1 : 0)
		ppfltData[nRow] = &ppfltData[nRow - 1][nWidth];

	// result
	return ppfltData;
}

void CScatBlockLP::freeData(std::complex<float>** pData, int nHeight, bool bSymetric)
{
	// number of complex we have
	int nSize = bSymetric ? (((nHeight + 1) * nHeight) / 2) : (nHeight * nHeight);

	// delete data by calling destructors
	while (--nSize >= 0)
		pData[0][nSize].~complex<float>();

	// delete all
	delete [] (BYTE*)pData;
}

/****************************************************************************
* CScatBlockLPT - scat matrix block transform for LP matrices				*
****************************************************************************/

/****************************************************************************
* Overrides																	*
****************************************************************************/
int CScatBlockLPT::getWidth() const
{
	// switch Transform and return width or height
	switch (m_nTransform)
	{
	case CTransform::eTransform_r180:
	case CTransform::eTransform_fr0:
	case CTransform::eTransform_fr180:
	return m_pOwner->getWidth();
	case CTransform::eTransform_r90:
	case CTransform::eTransform_r270:
	case CTransform::eTransform_fr90:
	case CTransform::eTransform_fr270:
	return m_pOwner->getHeight();
	default: // invalid
	assert(FALSE);
	}
	// invalid
	return -1;
}

int CScatBlockLPT::getHeight() const
{
	// switch Transform and return width or height
	switch (m_nTransform)
	{
	case CTransform::eTransform_r0:
	case CTransform::eTransform_r180:
	case CTransform::eTransform_fr0:
	case CTransform::eTransform_fr180:
	return m_pOwner->getHeight();
	case CTransform::eTransform_r90:
	case CTransform::eTransform_r270:
	case CTransform::eTransform_fr90:
	case CTransform::eTransform_fr270:
	return m_pOwner->getWidth();
	default: // invalid
	assert(FALSE);
	}
	// invalid
	return -1;
}

const CScatBlock* CScatBlockLPT::getTransformed(
	CTransform::eTransform nTransform, bool bCreateIfNotThere) const
{
	// validate params
	assert(nTransform >= CTransform::eTransform_r0 && nTransform < CTransform::eTransform_count);

	// now combine trasfroms to get the real one
	CTransform::eTransform nCombinedTransform =
		CTransform::combineTransforms(m_nTransform, nTransform);

	// return matrix reference from shared data
	return m_pOwner->getTransformed(nCombinedTransform, bCreateIfNotThere);
}

std::complex<float>& CScatBlockLPT::valueAt(int nRow, int nColumn)
{
	// switch Transform and return width or height
	switch (m_nTransform)
	{
	case CTransform::eTransform_r180:
	case CTransform::eTransform_fr0:
	case CTransform::eTransform_fr180:
	case CTransform::eTransform_r90:
	case CTransform::eTransform_r270:
	case CTransform::eTransform_fr90:
		// TODO: NOT FULLY IMPLEMENTED YET
		assert(FALSE);
	break;	
	case CTransform::eTransform_fr270:
	return m_pOwner->valueAt(nColumn, nRow);
	default: // invalid
	assert(FALSE);
	}
	// invalid
	return *(std::complex<float>*)NULL;
}

/****************************************************************************
* Construction/Destruction													*
****************************************************************************/
CScatBlockLPT::CScatBlockLPT(CScatBlockLP* pOwner, CTransform::eTransform nTransform) :
	m_pOwner(pOwner),
	m_nTransform(nTransform)
{
	// should never be null transformed
	assert(m_nTransform != CTransform::eTransform_r0);
}

/****************************************************************************
* End of file																*
****************************************************************************/
